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Abstract. We compute the Hugoniot curves of both neat TATB and its detonation products 
mixture using atomistic simulation tools. To compute the Hugoniot states, we adapted our "Sam- 
pling Constraints in Average" (SCA) method (Maillet et al., Applied Math. Research eXpress 2008, 
2009) to Monte-Carlo simulations. For neat TATB, we show that the potential proposed by Rai (Rai 
et al., J. Chem. Phys. 129, 2008) is not accurate enough to predict the Hugoniot curve and requires 
some optimization of its parameters. Concerning detonation products, thermodynamic properties at 
chemical equilibrium are computed using a specific RxMC method (Bourasseau et al., Phys. Chem. 
Chem. Phys. 13, 2011) taking into account the presence of carbon clusters in the fluid mixture. 
We show that this explicit description of the solid phase immersed in the fluid phase modifles the 
chemical equilibrium. 

I. INTRODUCTION 

To understand and model detonation phenomena, it is very important to describe accurately 
the thermodynamic properties of the neat explosive as well as its detonation product mixture. 
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When a material is hit by a shock wave, its thermodynamic state changes towards a state with 
higher temperature, higher density and higher pressure. The set of admissible states which can be 
attained from given initial conditions, for shocks of various strengths, is called the Hugoniot curve 
(see equation ([T]) for a precise definition). For explosives, exothermic chemical decompositions are 
triggered by the passage of a shock wave. However, the state of the material after the propagation 
of the shock is a particular point of the Hugoniot curve of the unreacted material, the so-called 
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ZND state, which corresponds to a high pressure and a relatively high temperature 
this state, exothermic chemical reactions produce small stable molecules (such as H2O, CO2, N2) 
called detonation products. As for every thermally activated processes, the initial thermodynamic 
conditions of the system before it reacts are one of the parameters controlling the reaction rate. 
Therefore, the determination of the temperature of the material after the passage of the shock 
wave, and more generally of the complete equation of state (EOS) of the neat explosive, is a key 
point for predicting the first step of the molecular decomposition. Moreover, the energy released 
by a given explosive can be estimated from the knowledge of the isentropic curve of the detonation 
products mixture in the pressure-volume plane. Thus, important efforts have also been dedicated 
in the past decades to measure and predict the equilibrium properties of the detonation products in 
the high pressure and high temperature regime. Unfortunately, measurements of thermodynamic 
properties of both shocked material and detonation product mixtures are relatively scarce. As a 
consequence, following a similar work on liquid nitromethane and other energetic materials % |5| , 
we compute here the Hugoniot curves of TATB and of its detonation products using techniques 
and models from statistical physics. 

To this end, we adapted the "Sampling Constraints in Average" method (SCA) proposed in |6j 
to our molecular Monte Carlo simulation tool. This method aims at sampling microscopic config- 
urations of a system, consistent with a thermodynamic ensemble (canonical, isobaric-isothermal, 
grand-canonical, etc), and such that some constraints are satisfied in average. A typical applica- 
tion is the determination of the temperature of a system given its average energy and its density. 
Here, we use SCA to compute the Hugoniot curves of both the neat explosive and the detonation 
product mixture. The composition at chemical equilibrium of the detonation product mixture. 
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n 

is obtained by the RxMC method 7]. It is important to take into account the presence of sohd 
carbon clusters in the detonation product mixture. The RxMC method models carbon clusters 
as mesoparticles embedded in the mixture and therefore allows us to calculate thermodynamic 
properties of a heterogeneous liquid/solid system at chemical equilibrium. 

This paper is organized as follows. In Section |TT1 we briefly review the method to sample con- 
straints in average, in the particular case of the computation of Hugoniot curves. We then describe 
the system we consider and the sampling method used to sample its microscopic configurations in 
Section IlIII Numerical results are then presented in Section llVl Our findings are finally summa- 
rized in the conclusion (Section IV)) . 



II. SAMPLING CONSTRAINTS IN AVERAGE 



A. Description of the method 



We present here the SCA method applied to the computation of Hugoniot curves when a Monte 
Carlo sampling method is used (see 6[ for a more general presentation and other possible applica- 
tions). 



1. Hugoniot curves 

The Hugoniot curve is the ensemble of accessible thermodynamic states that a system can reach 
from a given initial state after the passage of a shock wave. The thermodynamic quantities of 
a material in the initial unshocked state and the final shocked state are related by the Rankine- 
Hugoniot relation: 

E^Eo = ^{P + Po){Vo^V), (1) 

where P, V, and E are respectively the pressure, volume and energy of the system in the shocked 
state (at temperature T), and Pq, Vq, and Eq are the pressure, volume and energy of the system 
in the initial state (pole). 

Macroscopic relations such as ([T]) can be obtained by averaging functions of the microscopic state 



of the system, as predicted by the laws of statistical physics. In this study, we restrict ourselves 
to the NPT ensemble, but generalizations to other thermodynamic ensembles are straightforward. 
We denote by (g, V) the configuration of the system: q = {qi, . . . , g^) are the positions of the 
A particles, and V is the volume of the simulation box. The potential energy of the system is 
denoted by Epot{q)- A reformulation of ([l} consists in replacing E by the average microscopic 
energy obtained by averaging -Epot in the thermodynamic ensemble at hand and adding the kinetic 
contribution, and replacing V by the average volume at the given temperature and pressure. 
Therefore, ^ should be understood as 

{H)p,T - (i/)Po,To - l{P + Po){{V)p„To - {V)p,t), (2) 

where (•)p,t denotes an isobaric-isothermal average at fixed pressure P and temperature T. 
2. Previous methods to compute Hugoniot curves with Monte Carlo sampling method 



9|. In this method. 



Hugoniot curves were previously computed with the AE-EOS method 
the pressure P of the state on the Hugoniot curve is fixed, and the temperature such that ^ is 
satisfied is determined iteratively, using Newton's algorithm. More precisely, denote by T" the 
approximation of the temperature T at the n^^ step. A short NPT simulation at temperature T" 
and pressure P is run, and an approximation to the average 

i?g(T") = (if)p,T" -Eo-]^{P + Po) (Vo - (F)p,T") (3) 

is computed. The derivative of Hg with respect to T is then estimated at the temperature T" 
using a finite difference, and a new temperature is obtained as 

To initialize the algorithm, is chosen arbitrarily, and is obtained by adding a given AT to 
. The temperature T" converges to the temperature Th^, which corresponds to HgiTn^) = 0. 
The AE-EOS method gives good results, and allowed us to obtain the Hugoniot curves of various 
energetic materials and their detonation products mixtures. However, its major drawback is that 
the modifications of the temperature do not easily allow to perform a statistical ensemble average 
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on the full calculation. At the end of the AE-EOS calculation, a usual NVT or NPT simulation 
at temperature T — T^^ is needed to obtain the direct and derivative thermodynamic properties. 
Therefore, a lot of computational time is wasted in the process. 



3. Sampling constraints in average 

To get rid of this important drawback of the AE-EOS method, we have recently developed 
a more efficient method (called SCA) which allows to sample microscopic configurations with a 
constraint satisfied in average during the simulation. This method was developed in a molecular 

n 

dynamics framework ,6]. We extend it here to Monte Carlo simulations. Applied to the Hugoniot 
calculation case, it allows to determine the temperature Tu^ such that is indeed equal to zero in 
average, and it moreover allows to correctly explore the phase-space configurations consistent with 
the appropriate thermodynamic ensemble at this temperature. Thus, a single simulation gives 
the Hugoniot temperature and the thermodynamic properties of the system at this point. The 
convergence of this method was studied from a mathematical viewpoint in . 

Let us now present a brief description of the adaptation of SCA to a Monte Carlo sampling 
method such as standard Metropolis-Hastings algorithms 
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Heuristically, the idea is to 



select a control variable (the temperature T here), to decompose its values into bins, and to 
accumulate the instantaneous values of the property under investigation [H^ in our application) 
every time the control variable falls in the corresponding bin. This allows to construct estimates of 
the average property as a function of the control variable. To this end, a standard NPT simulation 
is performed, with an additional dynamical update of the temperature depending on the observed 
values of the average property which should be constrained. The key point is to determine how 
the temperature is changed. In our application, _ffg is positive when T > Th^ and negative when 
T < Th^- As our goal is to find the temperature T//^ for which is zero, it is natural to resort 
to the following temperature update: 

r"+i = ^"-a(i^g)p,T". 
In fact, the average value (iJg) p,t" is not known exactly, and should be replaced with its approxi- 



mation obtained from the histogram for the bin corresponding to T". It corresponds to an average 
computed along the whole trajectory, which is therefore more and more accurate as the simulation 
goes on. 

More precisely, consider a temperature grid = Tmin + iAT with i taken between and M , and 
rpM _ jj^^^ Denote by I(T) the function which returns the index of the bin corresponding to the 
temperature T, and by MC p_t{q, V) the Monte Carlo algorithm which returns a new configuration 
g', V' from a previous one, for a given pressure P and a given temperature T. The SCA algorithm 
then schematically reads as follows: 

n+l 



^ ( V") -Eo + hp + PoWo - V"')) Si^T^^^i^T 



T 



m=0 



n+l 

^/(T'")=/(T") 

rn—0 

See also Figure [T] for a cartoon representation of the procedure. In fact, for Monte Carlo sampling 
methods, it is convenient to repeat several times (say, N times) the update of the configuration 
before updating the temperature. 
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FIG. 1: Schematic representation of the SCA algorithm. 



It is expected that T" converges to . As for the continuous dynamics considered in 65 , such 
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a convergence can be expected only if the parameters of the algorithm are well chosen - although 
it is very easy to find satisfying values of these parameters, see Section IIIBI The two parameters 
that control the efficiency of the algorithm are a and N: 

(i) The quantity a controls the rate at which the temperature is updated. If a is too small, 
the temperature changes very slowly, leading to an inefficient algorithm. If a is too large, 
the variations are too brutal and the system is driven out of equilibrium. This may lead to 
numerical instabilities. A practical way of choosing a is therefore to start from a conservatively 
small value, and increasing it until some numerical instabilities are observed. This can be 
done using very short preliminary computations. Besides, our experience is that the precise 
value of a is not very important, and the range of admissible a is robust with respect to the 
changes in the thermodynamic conditions. In any case, the value given to a does not change 
the final result, it only has an infiuence on the efficiency of the convergence. 

(ii) The number of steps N between two temperature updates. Actually, in the original Molecular 

II 

Dynamics approach [6j], the temperature is changed at every time step. However, in a Monte 
Carlo simulation, the system may only be slightly modified at each Monte Carlo step, and it 
appears more convenient to let the system equilibrate a little while at a given temperature 
before changing it again. As a consequence, a reasonable value for N could be an entire 
Monte Carlo cycle between two temperature changes. 

B. Application to a model system 

To understand the infiuence of the parameters N and a, we performed several simulations 
with this method on a test system composed of 400 point-like particles interacting with a standard 
pairwise Lennard- Jones 6-12 interaction potential depending only on the relative distances between 
the particles. A cut-off equal to the half of the box length has been used to reduce the computing 
time. Periodic boundary conditions together with long range corrections have been used. 

Parameters have been taken arbitrarily equal to e = 120.0 K and cr = 3.40 A. We computed the 
temperature for which Hg = 0, at P = 1 GPa, with initial thermodynamic conditions arbitrarily 
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chosen as: Eq = 345 J.g"\ Vq = 5.848 10"'^ m^.g-i and Pq = 1 10^ Pa. Monte Carlo moves 
in these test simulations are translation moves (95 %) and isotropic volume changes (5 %). 
Translation distances and volume changes coefficients are selected randomly under maximum 
values that are updated every 10000 MC moves in order to obtain 40 % of accepted moves. A 
standard Metropolis scheme is used to accept or reject new configurations. Finally, around 1000 
MC moves are needed to complete a Monte Carlo cycle. 





Sampling Constraints in Average (New method) 


AE-EOS (Old method) 


N 


1000 


1000 


1000 


100 


10000 


10000 


100000 


a (K.g.J-i) 


1. 


0.1 


0.01 


0.1 


0.1 






Th, 


1285.7 


1284.0 


1284.9 


1284.7 


1284.8 


1279.2 


1285.6 


(K) 


± 0.4 


± 0.04 


± 0.06 


± 0.2 


± 0.2 


± 0.0 


± 1.8 


^^g 


-0.06 


0.13 


0.14 


-0.38 


0.06 


-2.22 


1.40 




± 0.53 


± 0.42 


± 0.57 


± 0.56 


± 0.61 


± 0.46 


± 1.53 



TABLE I: Comparison between results obtained on the test system with the SCA method (New method), 
and the AE-EOS method. 

Table U shows the results obtained for different values of a and A^. We also show in this table 
the results obtained with the AE-EOS method. In the two cases, the parameter N corresponds 
to the number of MC moves performed between two temperature changes. A total of 10^ AIC 
moves have been performed in all cases, and statistical errors have been estimated using blocks 
averages over the 5.10^ last iterations. Concerning the SCA method, it appears that the obtained 
temperature does not depend on the parameters, as expected. The "Hugoniot" temperature is 
equal to 1285 ± 1 K. Aloreover, the constraint iJg = is respected in average for all the parameter 
sets considered. Note that in such a system, the total energy is about several tens of J.g^^, so that 
the statistical uncertainties observed with both methods are very small. However, the results from 
the new method arc more accurate and reliable that those obtained with the AE-EOS method at 
a fixed computational cost. Moreover, it appears that AE-EOS method can converge to wrong 
values of the temperature: in the example corresponding to N = 10000, the temperature seems 
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extremely well converged with a very reduced statistical error, but with a corresponding Hg value 
not equal to zero. 




FIG. 2: Evolution of the temperature during the simulation for various values of the parameter a (in 
K.g.J"^). The value A'' is kept constant and equal to 1000. The insert corresponds to a zoom on the first 
lO" iterations. 

Figure [2] shows the evolution of the temperature observed during the simulation, for various 
values of a, when N is kept constant and equal to 1000. As shown in Table U the temperature 
obtained at the end of the simulation does not depend on the parameter a. However, too high a 
value of a triggers oscillations of the temperature around the target temperature, leading to higher 
statistical uncertainties. On the opposite, too small a value of a leads to a very slow and inefficient 
convergence. 

Figure |3] shows the evolution of the temperature observed during the simulation, for various 
values of N, keeping a constant and equal to 0.1 K.g.J^^. Once again, the limiting value of the 
temperature does not depend on the value of N. However, it seems that too high a frequency of 
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temperature updates (corresponding to a small value of N) can lead to instabilities, whereas too 
high a value of N slows down the convergence. 
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FIG. 3: Evolution of the temperature during the simulation for various values of the parameter A'^. The 
value a is kept constant and equal to 0.1 K.g.J^^. The insert corresponds to a zoom on the first 10® iter- 
ations. 



III. DESCRIPTION OF THE SYSTEM 



We present successively the Monte Carlo method used to sample the configurations of the system 
fSection llll Ap . and the atomistic description of the interactions f Section IIII B[) . 



A. RxMC with mesoparticle 

The goal of the RxMC method is to predict equilibrium properties of multi-component systems 
constrained by chemical equations. Microscopic configurations consistent with the chemical equi- 
librium can be sampled by resorting to a specific thermodynamic ensemble: the Reaction Ensemble. 
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RxMC is a well-established molecular-level simulation method to sample the Reaction Ensemble 
and it has found applications in a large variety of problems. Interested readers are referred to 
previously published articles on the method and its application. In particular, Turner et al. have 



published a review of the method and its applications prior to 2008 



more recent applications are J, [7 
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12[, and some references for 



14| . The complete definition of this ensemble, including an 



15j. J.K. Johnson has 



16| . Another 



expression of its probability density, was given first by Smith and Triska 
also given a very nice derivation of the Reaction Ensemble acceptance probability 
interesting early reference is the pioneering work of M.S. Shaw, who proposed a method similar in 
nature to RxMC to simulate chemical equilibrium of molecular mixtures 
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18|. 



The RxMC method is based on a particular Monte Carlo move, the reaction move. For a given 
chemical equilibrium, this move consists in changing the chemical composition of the system by 
deleting reactant molecules, and inserting products molecules in order to keep constant the number 
of atoms of each species. Of course, the move has to be performed in forward and backward 
directions to ensure micro-reversibility. The averaged chemical composition of the simulation 
corresponds to the equilibrium composition of the system. 

In the reference 7[, we give a detailed description of the original method, and we also describe 
the new method we have proposed to perform RxMC simulation of detonation product mixtures 
that include solid carbon clusters in the fluid phase. It is indeed very important to explicitly model 
this phase of carbon in order to obtain correct simulation results since the chemical equilibrium of 
a single heterogeneous system is thermodynamically different from considering two homogeneous 
separated systems {i.e. liquid and solid) in equilibrium. Significant differences have been found 
between thermochemical calculations (considering separated phases) and explicit microscopic sim- 
ulations of the heterogeneous system 7|. 

The method described in 7[ models the solid phase as a mesoparticle representative of a cluster 
of Nc carbon atoms (denoted by MPjVc in the sequel), immersed in the reacting fluid. In this 
case, the mesoparticle can be considered as a single rigid molecule, and the chemical equations we 
consider to take into account the chemical equilibrium between the fluid and the solid phases are 
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similar to the following example: 

2 CO + MP No ^ CO2 + MP jvc+i (4) 

In the forward direction, the corresponding reaction move consists in deleting two CO molecules, 
and the mesoparticle of Nc carbon atoms, and inserting a CO2 molecule and the new mesoparticle 
MPjVc+i- Of course, an efficient way to proceed is to replace the old mesoparticle by the new 
mesoparticle, which corresponds in practice to changing the mesoparticle volume according to the 
addition of a carbon atom. This move is no longer completely random, and the bias has to be 
corrected. This is easily done by performing deletions, insertions and mesoparticle replacements 
always in the same order (see [J] for further precisions). 

It is possible to write the expression of the acceptance probability of such a move from the 
standard equations of the RxMC ensemble (see ,7|, which is dedicated to the precise description 
of this method). Denoting by s the number of chemical species (where the sth species is the solid 
carbon), it holds 

Pace = mm I 1, {PolSVy ■ exp I I ' 

where Po is the standard pressure, /? ~ l/{kBT), V is the volume of the system, ^ is the move 
direction (^ = +1 if the reaction move is performed in the forward direction, or -1 if it is performed 
backward), i> — with i^i the stoichiometric coefficient of the species i involved in the 

reaction (for the solid carbon, the notation i/g = i^Csoi is used), A{G^{T) is the standard Gibbs free 
energy of formation of the zth species at temperature T, Af Gcg^i {T, P) is the Gibbs free energy of 
formation of a mole of solid carbon cluster at temperature T and pressure P, Ni is the number of 
molecules of the ith species in the system, and AJ7 is the energy difference between the old and 
new configurations. 

It is interesting to note that the only input data needed in the reaction move are the free energies 
of formation Af G° (T) of the different species (for i < s — 1) and Af Gcsoi {T, P) ■ AH these quantities 

tI- It is 



can be obtained from experimental databases or numerical studies 



is also possible to compute 



a chemical equilibrium involving several chemical reactions. To this end, we add a preliminary step 
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consisting of randomly choosing the chemical reaction before each reaction move. Finally, using 
in addition the usual volume change of standard MC simulations, it is also possible to simulate a 
chemical equilibrium at constant pressure. 

Combining translation, rotation and reaction moves, and volumes changes, it is possible to 
simulate a chemical equilibrium at a given temperature and a given pressure. Statistical biases 
can be resorted to in order to improve the acceptance probability of insertion in dense phases. 



We employed to the pre-insertion bias method 



19[, already implemented in our MC code. An 



alternative (and similar) technique is the cavity bias sampling method of [2d |. The biases applied 
to the reaction move consist in inserting the first product molecules at the empty locations of 
the previously deleted reactant molecules. If other product molecules have to be inserted {i.e. if 
D > 0), the insertion is performed at a pre-selected location. The details of the algorithm and the 
expression of the acceptance probability when using the pre-insertion bias are presented in 4|. 

B. Force fields 

To perform all atom simulations of neat TATB, two force fields are available in the literature, one 



from Gee et at, proposed in 2004 



23. The first 



and another from Rai et al., proposed in 2008 
force field is not suited for Monte Carlo simulations, since molecules are supposed to be flexible, 
which imposes some additional expensive internal relaxation moves. As intramolecular hydrogen 
bonds and benzene cycles impose a relatively high rigidity of the molecular structure, we believe 
that the rigid molecule approximation is valid in this case. Moreover, atomic charges given by [21 [ 



seem to be unreliable. The second force field describes rigid molecules and contains more realistic 
atomic charges. We therefore decided to use the potential given in [22] in our simulations. This 
potential has been developed from the TraPPE force field of aniline and nitrobenzene molecules, 
whose molecular structures are close to the TATB structure. Parameters of intermolecular atomic 
Lennard- Jones 6-12 potentials have been transferred directly from those two molecules whereas ab 
initio calculations have been used to set the molecular geometry and the atomic charges. Table HIl 
gives the complete list of parameters used to model the neat TATB molecule. 

The experimental molecular structure of TATB has been measured at ambient temperature by 
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a {A) 


6(K) 


q(e) 


C (-NO2) 


3.60 


30.7 


-0.242 


C (-NH2) 


3.60 


30.7 


+0.408 


N (-O2) 


2.90 


30.0 


+0.007 


N (-H2) 


3.25 


160.0 


-0.738 





2.70 


42.0 


-0.104 


H 


0.50 


12.0 


+0.386 



TABLE II: Lennard- Jones 6-12 parameters and atomic partial charges used to model the atomic force 
centres of neat TATB. 



Cady and Larson in 1965 23|. The experimental geometry is given in Table Hill Dihedral angles 
are equal to except CCNO which is equal to 12°. Experimentally, the TATB molecule is neither 
absolutely plane, nor absolutely symmetric. The parameters used in our simulations are given in 
Table ITVl All dihedral angles are set to 0. 



dx-Y 


(A) 


X Z (°) 






C-C(-NH2)-C 


117.9 


C-C 


1.442 


C-C(-N02)-C 


122.0 


C-N(-H2) 


1.314 


C-N-0 


121.0 


C-N(-02) 


1.419 


O-N-O 


117.9 


N-O 


1.243 


C-C-N(-02) 


119.0 






C-C-N(-H2) 


121.1 



TABLE III: Experimental molecular structure of TATB at ambient temperature 23 



We consider 8 different chemical species to represent the detonation products mixture of TATB 
(solid carbon, CO2, H2O, CO, N2, H2, NH3 and CH4). Molecules in the fluid phase have been 
modelled through the exp-6 potential from Fried et al. [24]. The solid carbon has been modelled 
through a mesoparticle representing a carbon cluster, using the model we developed recently 7|. 
In this model, the inner properties of the cluster are given by an equation of state. Radius of the 
mesoparticle, and interaction potential between the mesoparticle and a fluid particle have been 
obtained from molecular dynamic simulations of all atoms carbon clusters using the LCBOPII 
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dx-Y 




X Z (°) 






C-C(-NH2)-C 


118.91 


C-C 


1.437 


C-C(-N02)-C 


121.09 


C-N(-H2) 


1.317 


C-N-O 


120.66 


C-N(-02) 


1.422 


0-N-O 


118.68 


N-O 


1.235 


C-C-N(-02) 


119.45 


N-H 


1.014 


C-C-N(-H2) 


120.54 






H-N-C 


116.53 






H-N-H 


126.94 



TABLE IV: Molecular structure of TATB used in our simulations 



potential 



251 . Details are given in reference 



7|. We also show in this reference that this model is 



particularly well suited to represent the effect of a real carbon cluster immersed in a fluid mixture. 



IV. NUMERICAL RESULTS 



The preliminary study presented in Section FlI Bl gives an order of magnitude for the parameters a 
and N . We performed a few small SCA preliminary runs for the TATB model at hand, and decided 
to use the following parameters to compute the Hugoniot curves: for neat TATB, N ~ 2000 and 
a = 1. K.g.J^^, while for the detonation product mixture of pure TATB, N = 20,000 and a — 
0.2 K.g.J-i. 



A. Hugoniot curve of neat TATB 



We first computed the thermodynamic properties of neat TATB at the pole, which corresponds 
to To = 300 K and Pq — 10^ Pa. The results are given in Table fVl They show that the computed 
direct thermodynamic properties and cell parameters are close to both experimental results and 
Rai's calculations. Differences between our results and Rai's results mainly arise from the way 
the long range electrostatic interactions are taken into account: we used the reaction field method 
28j whereas 



22l | used the Ewald summation method. On the other hand, there are important 
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Experimental 


Rai [22] et at 


This work 






Results [23] 


[22] 






Po (g-cm ^ 


1.938 


1.929 


1.932 


Direct 


Po (Pa) 


1.10^ 


1.10^ 


1.10^ 


Properties 


Vo (m^g"i) 


5.16.10"^ 


5.184.10"^ 


5.175.10-'^ 




Eo (J.g-') 






1346.7 




a (A) 


9.010 


9.05 


9.037 




b(A) 


9.028 


9.04 


9.018 


cell 


c(A) 


6.812 


6.80 


6.802 


parameters 


a(°) 


108.59 


110.0 


110.087 




HI 


91.82 


88.2 


88.086 




in 


119.97 


120.2 


120.185 




Cp (J.K-\g-i) 


1.00 [26] 


- 


0.665 ± 0.01 




Kg (GPa) 


17.3 [26] 




9.907 ± 0.29 


Derivative 




18.9 [27] 






Properties 


r 


0.20 [26] 




1.357 ± 0.12 




Cs (m.s-^) 


1460 [23] 




2338.8 ± 44.6 






304 [26] 




163.03 ± 7.6 



TABLE V: Direct and derivative thermodynamic properties and cell parameters of TATB calculated at 
ambient temperature (pole properties), and compared to Rai's results 23] and various experiments. Cp: 
calorific capacity at constant pressure, Kq: bulk modulus, F: Gruneisen coefficient, Cs: sound velocity, 
av- thermal expansion coefficient. 



differences in the derivative properties between the computed values and the experimental results. 
Although the reliability of experimental measurements is not completely guaranteed, the magnitude 
of the differences underlines a weakness of the potential, and shows that the potential does not 
reproduce correctly the evolution of the pressure as a function of the volume (Kq), of the energy 
as a function of the temperature (Cp), and of the volume as a function of the temperature {ay) 
near the initial state conditions. 

To analyze the behavior of the potential under pressure, we performed 4 NPT simulations at 



300 K at different pressures P (2, 5, 10 and 15 
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GPa). Results are displayed in Table IVll and 



26| and Stevens 



27| in Figure HI The quoted 



compared to experimental results from Olinger 
experimental results have been obtained by compression of TATB powder, which roughly corre- 
sponds to hydrostatic compression of the monocrystal. According to Stevens, discrepancies ob- 
served between the results from the two authors come from the hypothesis made to obtain the cell 
parameters. Figure 2] shows that numerical results are in good agreement with experimental ones. 
Rai et at have tested their potential up to 7 GPa. We show that this potential gives satisfying 
qualitative results up to 15 GPa. A careful investigation however reveals that the curvature of the 
P-V isotherm is not well reproduced, the curve being too convex. This is consistent with the fact 
that the computed bulk modulus Kq is smaller than the experimental one. 



P (GPa) 


V/Vo 


V (cm3.g-i) 


10-* 


1.0 


0.517 


2.0 


0.901 


0.467 


5.0 


0.843 


0.436 


10.0 


0.793 


0.410 


15.0 


0.764 


0.395 



TABLE VI: Calculated volumes along the isotherm 300 K. Statistical uncertainties are under 0.1 %. 



To calculate the Hugoniot curve of neat TATB, we performed 4 simulations using the SCA 
method described in Section at P = 2, 5, 10 and 15 GPa. Results are shown in Tabl e I VIII and 

301. The 



Figure [5l and are compared to experimental results from Coleburn et al. [29] and Marsh 
latter Hugoniot experimental measurements were obtained from compressed powders slightly less 
dense (1847 kg.m^'^ and 1876 kg.m^^^) than monocrystals (1938 kg.m^^^). Figure [S] shows that the 
pressure along the Hugoniot curve for a given compression is slightly overestimated. The difference 
seems sufficiently small to be corrected by a modification of the parameters of the potential using 
some potential optimization, as we previously did for nitromethane 3l|. Table IVlIl demonstrates 
the accuracy of the SCA method for obtaining Hugoniot states. The statistical uncertainty on the 
computed temperatures is really small, and the average value of Hg is close to zero, knowing that 
the total energy of such system is typically always over 1000 J.g~^. 
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FIG. 4: TATB isotherm 300 K: pressure as a function of specific volume. Calculation results are compared 
to experimental ones from Ofinger 2Q| and Stevens 27 1. 



P (GPa) 


V/Vo 


V (cm^g-i) 


T(K) 


Hg (J.g-') 


10-* 


1.0 


0.5175 


300.0 




2.0 


0.903 


0.4675 


342.9 ± 0.14 


-0.048 ± 0.15 


5.0 


0.8458 


0.4377 


403.9 ± 0.11 


0.068 ± 0.22 


10.0 


0.7962 


0.4120 


520.3 ± 0.11 


-0.001 ± 0.25 


15.0 


0.7659 


0.3963 


652.9 ± 0.17 


-0.045 ± 0.36 



TABLE VII: Pressure, compression, volume and temperature of TATB along the Hugoniot curve. The 
last column gives the average value of Hg calculated during the simulation. Statistical uncertainties on 
volumes are under 0.1 %. 
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FIG. 5: Hugoniot of neat TATB. Calculation results are compared to experimental ones from Coleburn 
and Marsh S^l- 

B. Hugoniot curve of detonation products of TATB 

The detonation product mixture of TATB is mainly composed of 8 molecular species (solid 
carbon, CO2, H2O, CO, N2, H2, NH3 and CH4). To model the global chemical equilibrium 
occurring in the system, the 4 following independent elementary chemical equilibriums have been 
considered simultaneously. Those 4 chemical equations have been determined using the Smith and 
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Missen method '33] , explained in 



2 CO 
2NH3 
CO + 2NH3 
CO2 + H2 



CO2 + Csolid 

N2 + 3H2 

N2 + CH4 + H20 

CO + H2O 



We computed 6 Hugoniot states of the system at P = 15, 20, 25, 30, 35 and 40 GPa, using 
the SCA method in the RxMC ensemble. The initial conditions {Eo,Po, Vq) are the same as for 
the unreacted case. These computations were performed under two different hypotheses: either 
considering the two phases completely separated (using RxMC with the Composite Ensemble 9|), 
or modelling the solid^hase as a mesoparticle immersed in the fluid phase (using the RxMC with 
mesoparticle method [7]). Results are presented in Table IVlIIl and in Figure [S] 





Separated phases 


Mesoparticle 


P (GPa) 


V/Vo 


V (cm^g-^) 


T (K) 


Hg (J.g-^) 


V/Vo 


V (cm^g"') 


T (K) 


Hg (J.g-^) 


15.0 


0.9443 


0.5023 


2821.5 ± 0.2 


0.358 ± 3.5 


0.9739 


0.5180 


2596.4 ± 4.9 


-0.035 ± 7.9 


20.0 


0.8703 


0.4629 


2902.8 ± 0.4 


-0.350 ± 3.8 


0.8942 


0.4757 


2657.0 ± 1.4 


-0.517 ± 5.7 


25.0 


0.8202 


0.4363 


2992.7 ± 0.4 


0.775 ± 4.0 


0.8408 


0.4472 


2720.1 ± 3.2 


-0.268 ± 3.6 


30.0 


0.7829 


0.4164 


3093.4 ± 1.9 


-0.893 ± 4.9 


0.8017 


0.4264 


2785.5 ± 5.9 


0.997 ± 5.3 


35.0 


0.7536 


0.4008 


3208.5 ± 1.4 


-0.857 ± 4.7 


0.7706 


0.4099 


2867.7 ± 5.6 


0.199 ± 4.4 


40.0 


0.7294 


0.3880 


3339.1 ± 4.3 


-0.771 ± 4.8 


0.7466 


0.3971 


2936.0 ± 3.9 


0.399 ± 3.8 



TABLE VIII: Pressure, compression, volume and temperature of detonation product mixture of TATB 
along the Hugoniot curve in two different cases: either the two phases are considered as completely sepa- 
rated, or the solid phase is modelled through a mesoparticle immersed in the fluid phase. The last column 
gives the average value of Hg (see equation [S]) calculated during the simulation. Statistical uncertainties 
on volumes are under 0.1 %. 



Table I Villi shows that the SCA method gives also good results in this case: The values of Hg 
are really close to zero. The uncertainties on the temperatures are around a few Kelvins, and 
uncertainties on Hg are around a few J.g~^, what appears satisfying knowing that the energy of 
such systems is typically superior to several hundreds of J.g^^. This is clearly satisfying, even if 
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• Mesoparticle 
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FIG. 6: Hugoniot of detonation product mixture of TATB in two different cases: either the two phases 
are considered as completely separated, or the solid phase is modelled through a mesoparticle immersed 
in the fluid phase. 

uncertainties are higher than for neat TATB. This is due to the fact that there are intrinsically 
more fluctuations on the detonation product mixture than for neat TATB since the system is in a 
fluid phase, and its chemical composition varies. 

Figure [6] shows that significant differences (up to 4.3 % on the calculated volume at a given 
pressure P) appear between the results obtained with the two different assumptions on the mod- 
elling of the solid phase (separated phases vs. mesoparticle immersed in the fluid phase) . This has 
already been shown and explained in a previous paper [tI . The discrepancies can be attributed to 
the expansion of the solid phase, which needs more energy in the second case because mesopar- 
ticles have to overcome the fiuid pressure to grow. As a consequence, the chemical equilibrium 
of the system is displaced towards less carbon atoms in the solid phase for a given pressure and 
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temperature. The consequence of this difference on the Hugoniot curve is difficuh to explain, but 
Figures [7] and |S] show that the evolution of the temperature and the amount of carbon in the solid 
phase along the Hugoniot are really different in the two situations. The differences in the computed 
temperatures with the two hypotheses can reach 12 % at high pressures. The difference in the 
amount of carbon atoms in the solid phase can reach 30 % at 40 GPa, and we also observe that 
the qualitative evolution of the amount of solid carbon is totally different. This is a supplementary 
evidence that the heterogeneity of the system must be taken into account. 
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FIG. 7: Pressure/Temperature evolution along the Hugoniot curve of detonation products of TATB in two 
different cases: either the two phases are considered as completely separated, or the solid phase is modelled 
through a mesoparticle immersed in the fluid phase. 
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FIG. 8: Evolution of the percentage of carbon atoms in the sohd phase along the Hugoniot of detonation 
product mixture of TATB in two different cases: either the two phases are considered as completely 
separated, or the solid phase is modelled through a mesoparticle immersed in the fluid phase. 

V. CONCLUSIONS 

This work has shown that the "Sampling Constraints in Average" method can be used to cal- 
culate, in the same framework, the Hugoniot curves of neat TATB and of detonation products 
of TATB, with efhciency and reliability. The statistical uncertainties obtained both on converged 
temperatures and average values of Hg are more than reasonable. Besides, the parameters of the 
method are easy to determine. Moreover, derivative thermodynamic properties can be computed 
accurately. 



The potential we used to model the neat TATB, proposed by Rai et al. 231, allows to reproduce 
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correctly the direct thermodynamic properties of TATB at the pole conditions. Nevertheless it fails 
to reproduce quantitatively the derivative properties, and in particular the compressibility of the 
system. As a consequence, the 300 K isotherm and the Hugoniot curve are both too convex. This 
could probably be corrected with an appropriate modification of the parameters of the potential 
through some optimization procedure 3l|. This work is currently in progress. 

We used the specific RxMC method 7] to obtain the chemical equilibrium of the detonation 
product mixture of TATB, including solid carbon clusters. In this method, the all-atoms carbon 
clusters are replaced by mesoparticles. Combining this version of the RxMC method and the SCA 
method allowed us to compute the Hugoniot curve of the detonation product mixture of TATB. 
As expected, the results differ significantly from the results obtained with the composite ensemble, 
where the two phases are considered as completely separated. This is an evidence of the fact 
that the heterogeneity of the system (i.e. the fact that the carbon clusters are immersed in the 
detonation products fluid) must be taken into account. This is particularly important to calculate 
the amount of carbon atom included in the solid phase. 

To conclude, let us emphasize that the detonation velocities computed from the results of this 
work cannot be compared with experimental detonation velocities because the equation of state 
used to model the solid phase is not representative of a cluster phase. Nevertheless, we are currently 
working on the parametrization of a carbon cluster equation of state based on molecular dynamics 



simulation results obtained with the LCBOPII potential 



25|. Using this new equation of state, we 



anticipate a more accurate prediction of detonation velocity of TATB. 
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